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1. INTRODUCTION 

Modeling of the electrical circuits which contain passive elements 

(R, L, C) and active nonlinear elements (transistors, diodes) , 

allows for only very simple circuits to be solved using analytical 

s t ncl 

procedures. To solve an electrical network Kirchhoff's I s or 2 n 
law (node analysis or mesh analysis, respectively) and Ohm's law 
are used. More than three nodes or meshes in the linear case 
requires a complex analytical solution since it leads to a system 
of algebraic equations with three or more unknowns. This fact 
restricts the use of the analytical method. Further, the nonlinear 
elements cannot be included in the network, because even the most 
simple circuit with nonlinear elements - a series connection between 
the resistor and the diode - cannot be solved by the analytical 
procedure. Therefore, computers using numerical analysis procedures 
are needed. This acceptance of the computer as a tool allows for 
an extreme increase in both the possibilities and demands. The 
size of the network may be up to or exceeding 500 nodes or meshes, 
and the nonlinear elements are included with all their nonlinear 
characteristics . 

Assessment begins with a description of the network, an understanding 
of which elements it may contain, and the selection of the proper 
mathematical procedure to use. Elements which the network can 
contain are resistors, inductors, capacitors, transistors, (BJT, 

JFET, MOSFET) , diodes and voltage and current sources (constant, 
time varying and those which depend upon some other branch voltage 
or current) . 

Thus, both linear and nonlinear elements are included. Linear 
elements are the constans, and the nonlinear are described as 
I = f(V) of V = f(I). An example of a nonlinear element is the 
diode, where the diode current is equal to 

2V D 

I D " I s (e kT “ 1) 
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The transistor cannot be described in such a simple way (with one 
equation) . This in done with a number of equations which lead to 
a transistor model. 

The nonlinear elements are therefore described by their model. 

For example, the bipolar transistor can be interpreted with Ebers 
Moll's model as shown in Fig 1. 



I t = I, s (e ^-l) 

I. = I.s (e^-l) 
Coe= Z N dC„I e -£r 
C D £ rTxoCr Ic T?T 

r ^tco 

T ‘ CT 

V 1 V* 


Crt = 


\V~ 


V* 

^ TEO 

Vbe 
V* 


Fig. 1. Ebers-Moll* s model of the 
bipolar transistor . 

After the network elements have been described the voltage and 
current sources are given. The sources can be fixed, time varying 
and dependent. For example, in transistors, a collector current 
generator is depedent upon the emitter current so, 

I C = a N I E ' 

and the generator I c is a dependent current source. 

After the electrical network is described, a choice of the 
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necessary calculation is made. The following types of analysis 
can be done : 


1. An operating point of the circuit, so called DC analysis - 
analysis of the steady state conditions. 

2. A response to the small alternating signal, so called AC 
analysis which is executed so that first the DC solution of the 
network (we calculate the operating point for all nonlinear 
elements of the network) is found i.e. we linearize every nonli 
near element about its operating point. 

The linearization of a nonlinear element is shown using a diode 
as an example (Fig. 2) . 

I„ 


Fig . 2. Diode characteristics . 

Let V. , Ij be the diode operating point. The linearization is 
obtained so that, instead of the nonlinear characteristics 



I 


D 



1 ) 


the tangent in the point V^, 1^ is used : 



V = V 


1 


which gives the conductance 



d I n 

G = ■ ■■ = H_ i 
^ d V D kT "-S 


e <3 1 = L I 
e kT kT D 


V = V, 


so, in an electrical network the diode is presented only with 
the conductance G. 


After linearizing the network, the small alternating signal 
is applied and the response at a given frequency f is calculated. 
In this type of analysis the response inside the frequency band- 
width from f i to f 2 is of interest. In this case, the network 
calculation is repeated several times, and the distance of the 
frequencies between f.^, f 2 is arbitrarily divided. Usualy, the 
geometrical division is taken, i.e. 

f = q m f 1 (f x , qf 1# q 2 f 1 ... q^ = f 2 > . 


For example, to calculate the response of the video amplifier 
between 1 Hz and 1,05 MHZ with q=2 gives : 


f ■ 1, 2, 4, 8, ... , 524288, 1,048576. 

In the given example, the calculation will be repeated 20 times. 

As a result, the frequency and phase amplifier characteristics 
are obtained. 

3. The third type of analyisis is the time domain analysis or 
transient analysis. In this type of analysis the DC network state 
is calculated and then the analysis is begun. This type is the 
most complex and involves the use of elaborate mathemetical and 
programming procedures. 

Nonlinear elements are presented with their linear models as in 
Fig 3. 

After the nonlinear elements are linearized, the network, which now 
represents a number of resistors and fixed voltage and current 
sources, is solved. When the DC condition is determined, due to 
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the reactive elements C,L and the changeable signal source, the 
state of the network must be forseen at some time t = t + At. 
After this is accomplished (this will be covered later) , the 
network is linearized again, but now with the new conditions, 
since voltages and currents through the nonlinear elements have 
changed. This procedure is repeated step by step until time 
t = t gtOD is reached, at which time the analysis ends. 



Fig. 3. Linearized diode model. 
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2. NETWORK ANALYSIS NODAL APPROACH 


st 

Network analysis can be performed using either Kirchhoff's 1 
or 2 Law (nodal analysis or mesh analysis, respectively). The 
nodal analysis is more convenient and will be applied below using 
a simple resistive network (Fig. 4). 



Fig . 4. Oriented linear network. 

The network description must be done in a simple and definite way. 
Nodes and branches of the network are assigned symbols. One node 
is called the reference node 0. In the given example the network 
is described as in Fig. 5. 


R1 101k 



IS1 0 1 1 A JI-Ij 

IS5 3 2 .1 A v 


Fig. 5. Network description and 

network branch definition. 
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The given description in Fig. 5 completely defines the network 
so that the network may be reconstructed and solved. Solution of 
the various equations is done by computer analysis. Obviously, 
such a program is large and has 2000 - 10 000 instructions. 


2.1. Conductance Matrix 


After the network is described all the relevant equations are 
written. All mathematical equations should be in the matrix form, 
because this form is very convenient for mathematical manipulation 
and programming. 

For the network in Fig. 4, the following nodal equations are 
written : 


node 

1 

i 1 +o+o+i 4 +o = 

X sl 



node 

2 

o + i 2 + o - i 4 + i 5 = 

I s5 


(1) 

node 

3 

0+0+I 3 +0-I 5 = 

-I s5 



further 






node 

1 

G 1 V 1 + G 4 (v i " 

V 

= 

X sl 

node 

2 

” G 4 ( V i- V 2> + G 2 V 2 +C 5 

(V 2" V 3 ) 

= 

T s5 

node 

3 

-G 5 <v 2 -v 3 ) + G 3 V 3 


= 

_I s5 

rearranged 



\ 


node 

1 

(G 1+ G 4 )V 1 - g 4 V 2 + 0 


= 

X sl 

node 

2 

— G 4 v 1+ (g 2+ g 4 +g 5 )v 2 -g 5 v 3 

= 

I s5 

node 

3 

0 - G S V 2 +(G 3 +G 5 )V 3 


= 

_I s5 


( 2 ) 
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In the matrix form, term (2) is given as 


node 1 

G l +G 4 ~ G 4 0 


V 1 



node 2 

-G 4 G 2 +G 4 +G 5 ~ G 5 


V 2 

= 

X s5 

node 3 

0 " G 5 G 3 +G 5 


V 3 


“ I s5 


node 1 node 2 node 3 


( 3 ) 


or, in a shortened form : G V = I 

9 ^ — n — sn 

where G is a conductance matrix from (3) , vector V is the 

r \/ ””11 

voltage vector of the nodes towards the reference node (node 0 ) , 
and I on the right side is the current vector. Matrix G is made 
from the elemental network (Fig. 4) as follows : 

- a diagonal element ii is the sum of all the 
conductances that are connected with node i, 

- other elements in the rov/ are on the places •£ j , so 
that on the place ij comes the conductance which 
connects node i with node j , but with a negative 
symbol , 

- conductances, that are from node i connected, 
with a referent node, are only in the diagonal. 

The current sources vector element represents the sum of all 
current sources which enter that node. This procedure is very 
convenient for arranging the matrices from the topological 
network descriptopn which is described in Fig 5. Further, the 
relation G V = I mav be derived using Kirchhoff's nodal 
analysis, Ohra's Lav/, and the definition of a voltage drop on a 
resistor as demonstrated below. 


1 Kirchhoff's Lav/ . A standard network branch is defined 
in Fig. 6. Between nodes rn and n is a conductance G ; 
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parallel to it is a current source, I . 

s 



Fig. 6. Standard network branch. 


The total branch current is J = I - I .A network branch can 

s 

contain only the condudtance G, the current source I s , or both 
at the same time. Considering this definition, term (1) can be 
written as 


or as 


1 

0 

0 


0 

1 

0 


0 1 

0 -1 

1 0 



0 


1 

0 

0 


0 0 10 

10-11 
0 10-1 



or, in the matrix form 

A J = 0 


= 0 
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A (I - I 3 ) = 0 

A I = A I , 1 st It.L. (4) 

~ — — s 

Matrix A iz called the incidence matrix and can be constructed 

A/ 

iron a topological network description (Fig. 4) . Progressing 
from branch 1 to branch 5, a value of 1 is written in nlace of 
the node where the branch current exits, and -1 where the 
current enters. 

Matrix A is the connlete incidence matrix, and also contains the 
reference node: 


branch 

1 

2 

3 

4 

5 

node 

1 

1 

0 

0 

1 

0 

node 

2 

0 

1 

0 

-1 

1 

node 

3 

0 

0 

1 

0 

-1 

node 

4=0 

-1 

-1 

-1 

0 

0 


The fourth rov; of matrix A raav be reconstructed from the first 

~a 

three and is therefore unnecessary. For example, in the first 
column (branch 1) we find only +1. This means that -1 is in the 
fourth row, signifying that the first branch is connected between 
node 1 and the referente node 0. 


Fig. 4 indicates that the network has four nodes. Since the zero 
value may be associated to any single node, the real number of 
unknowns is three. 

Ohm's Lav; . For every network branch the connection between 
the voltage and the current may be written as follows : 
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G 1 v bl = *1 

G 1 


V bl 


*1 

G 2 V b2 = *2 

G 2 0 


V b2 


Z 2 

G 3 V b3 = *3 => 

G 3 


V b3 

= 

I 3 

G 4 V b4 = *4 

0 G. 

4 


V b4 


Z 4 

G 5 V b5 = *5 

G 5 


V b5 


X 5 


G d V b = I , Ohm's law. 


where G D is the diagonal conductance matrix, vector is the 
branch voltage vector, and I is the branch current vector. 


Branch Voltage Drop Definition . The voltage drop across 
a resistor is obtained by calculating the branch voltage from the 
difference between the potentials of the branch ends: 


V 1 * V bl 

10 0 


V 1 


v bi 

IT 

CM 

A 

> 

II 

CM 

> 

0 10 


V 2 

r= 

V b2 

V 3 = V b3 

0 0 1 


V 3 


V b3 

V 1“ V 2 = V b4 

1-10 




V b4 

V 2“ V 3 = V b5 

0 1-1 




V b5 


or 


£ 2 n - s. • 


( 6 ) 


A matrix which connects node and branch voltages is the 
transposed incidence matrix A^. 

Now we connect the terms (4), (5), and (6) 
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ii mi in mu ■ 


I !■ III I II 


O 


b I - h Is 

Sd — b ~ — r 


A V = V, . 
rJ — n — b 


A simple insertion gives: 


(7) 


A G_. V, = A I , 
r- — b — s ' 

A G n A t V = A I f 


(8) 


where G = A A fc . The conductance matrix is obtained using 
the incidence matrix and the diagonal conductance matrix. 

From term (8) , node voltages are calculated using : 





A I 
i s 


(9) 


All other parameters can be calculated from the known node 
voltages using term (7) : 

branch voltages = A*" V n , 

branch currents I = G_ V, = G^ A fc V 

— •~'D — b ~D ^ — n 

The matrix G = A G^ A fc for the network in Fig. 4 is calculated as: 

.-s' n-r rw»]j /*/ 


A G A 

-V rv 


1 0 
0 1 
0 0 


0 1 


0 


G 


1 


10 0 


0-1 1 


0 


0 10 


1 0-1 


0 


G 


3 


0 0 


1 


G 


4 


1-1 0 


G 


5 


0 


1 -1 
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G 1 

0 

0 G 4 

0 


1 

0 

0 


<G]+G 4 ) 

- G 4 . 0 

0 

G 2 

0 -G 4 

G 5 


0 

1 

0 


- G 4 

(G 2 + G 4+ G 5 ) - G 5 

0 

0 

G 3 0 

” G 5 


0 

0 

1 

= 

0 

(g 3 + g 5 ) 






1 

-1 

0 







0 

1 

-1 




2 . 2 Dependent Current Sources 

Any network branch can contain a current source. Current sources 
may be either dependent or independent. Independent current 
sources are placed on the right side vector as already shown. 
Dependent current sources are treated as follows : the source 
is described using the independent (control) branch voltage 
and is placed on the right side vector. Next, the terms are 
switched from the right to the left side (Fig. 7) . 

The independent branch current is given as 


h - G 1 (V ml - V m2' • 


and the dependent current source is 


I 2 = 6 I, - e Gi <V ml - v m2 ) 


This kind of dependent current source is called a current 
controlled dependent current source. If the source dependens 
upon the voltage, it is called a voltage controlled dependent 
current source and is expressed as 


- G m (v ' 


ml 


- v r«2> 
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independent branch 
(control branch) 

dependent branch 


independent branch 


current sources vector 




Fig. 7. Inserting of the dependent current 
source into matrix G. 

The dependent current source is determined by four terms in the 
matrix G which are not symmetrical with regard to a diagonal. Thus, 
the matrix G is generally symmetrical; however if there are some 
dependent current sources, matrix G, will have some nonsymmetrical 
terms. In this case, dependent current sources represented by 
following matrix G terms: 
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o 

c 

to 

u 

eQ 

-P 

c 

(1) 

T3 

C 

0) 

Pi 

Cl) 

TO 



BG 1 

- B G X 

or 

G 

m 

- G m 


- BG 1 

t 

BG 1 

t 


_G m 

t 

G m 

t 


control branch control branch 


2.3. Voltage Sources in the Network 

A voltage source which is alone in the branch / cannot be directly 
included in the network which is analysed by the nodal method. 

This is not possible because, by Ohm's law, term (5) , makes it 
impossible to calculate the current of the branch. 

Observing a network which has a voltage source, it is obvious 
that the number of the unknown node voltages is decreased by one. 
The known voltage source V , if attached between nodes V, 1 and 

V k2 ' gives V kl = V k2 + V ° r V k2 = V kl " V This rec 3 uires a 

smaller number of unknowns in the case of the network which is 
being solved by the nodal method. Voltage sources in the network 
are of two basic types: sources connected with one pole to the 
referent node and sources connected between any two nodes. 

a) Voltage sources connected with one pole to the 
referent node 

Most voltage sources that are applied in the network are grounded 
with one pole. These are supply voltages and input signal voltages 
First described is the insertion of these sources in the matrix 
term (8) . We will start from the example in Fig. 8. The network 
is described below assuming that the voltage source current (I x ) 
is known : 


15 



2 



Fig. 8. A network with a voltage source with 
one pole grounded. 


G l +G 2 

~ G 1 

0 

-G 2 


V 

s 


I 

X 

- G 1 

G l fG 3 

“ G 3 

0 


V 2 


0 

0 

“ G 3 

G 3 +G 4 +G 5 

” G 5 


V 3 


0 

- G 2 

0 

- G 5 

G 2+ g 5 


V 4 


0 


where V, = V . 

l s 

Switching the first column to the right the arrangement of the 

unknowns V n is changed from V g , V^, to V^, V^, V g , i.e. 

the known voltage V is moved to the bottom of the vector V . 

s -n 

After this operation, the matrix £ is nonsymmetrical and the 
current I has remained in the same place. The second operation 
switches the entire first equation to the last place, so that 
the arrangement of the current vector I is changed: I x , 0, 0, 0 
becomes 0, 0, 0, I . Aftr this, the matrix has become symmetrical 

X 

again and the unknown current I is placed at the bottom of the 
current vector. When both operations are carried out, term (10) 
becomes: 
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g 1+ g 3 


-G. 



" G 3 0 ! " G ! 

| 


V 2 

G 3 +G 4 +G 5 -°5 ! 0 


V 3 

-g 5 g 2 +g 5 , -g 2 

1 


I 

1 

> 1 
1 
1 

1 

0 G_ I G.+G„ 


V 

2 12 


s 


0 

0 

0 


(ID 


I 


x 


This term can be written in the condensed form as 


G 

G 


V 


I 

~nn 1 
J 

^ns 


— n 


-s 

l 

G i 

G 


V 


I 

vsn j 

r~SS 


— s 


—X 


From term (12) we get: 


( 12 ) 


G V + G V = I 
vnn — n ~ns — s — s 


G V + G V = I 
vsn — n v ss — s — x 


The first equation from (13) gives: 


G V„ = (I 
r^nn — n — -s 


G 

^ns 




and further 


V = G -1 (I 
— n <"nn — s 


G V ) . 
~ns —s' 


(13) 


(14) 


Thus the described procedure eliminates a number of the unknov/ns 
in the system and includes voltage sources. The second equation 
from (13 ) , after calculating the node voltages V n from (14) , 
enables us to calculate currents I which flovr through the 
voltage sources. 
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3y rearrangement of term (11) using (13) and (14) we obtain: 


G l +G 3 

- G 3. 

0 


V 2 


G. V 
1 s 

- C 3 

G 3 +G 4 +G 5 

"* G 5 


V 3 

= 

0 

0 

" G 5 

G 2 +G 5 


V 4 


G 0 V 
2 s 


k = -°1 V 2 + G 2 V 4 + <G 1 +G 2> v s 


If the network contains other voltage sources with one grounded 
pole, they are eliminated by switching all columns which contain 
the voltages V and the rows which contain the relevant currents 
I r to the right side, i.e. to the bottom of the matrix term. 

X 


b) Voltage sources connected between any 
two nodes 

In this case, the procedure is more complicated and is rarely 
applied in the electrical netv?or;- nodal analysis. Because of this 
the entire algebraic procedure will not be given. The basic prin- 
ciple involved is covered briefly. For a voltage source such as 

that in Fig. 9, the network, where the voltage source V is place 

s 

can be described as if it consisted of two separate voltage sour- 
ces V . After this, there are two voltage source branches and 
node is eliminated. The series connection between the voltage 
and the conductance is changed regarding Norton's theorem with 
a parallel connection between the current source and the conduc- 
tance as shown in Fig. 9. 


IS 







3. NUMERICAL SOLUTION OF ALGEBRAIC EQUATIONS 


3.1. Gaussian Elimination 


Chapter 2 shows how matrix G and vector 1^ are obtained, thus 

the network is described bv G V - I . The node voltages 

* — n -sn 

solution is symbolically v/ritten as 


V 


— n 





— sn 


(15) 


where G _1 is an inverse matrix of G. The solution will be explai- 
ned as it is the custom in mathematics in the example 


A x = b 


or, rewritten as 


x = A, -1 b . (16) 

The solution by Cramer's rule requires n! operations and cannot 

1 Q 

be used for larger matrices. (10! = 3,628,800, and 201=2.432*10 ) . 

Because of this, the system is solved by the Gaussian elimination 
method, or the modified Gaussian method which is sometimes called 
Crout's reduction. This procedure requires about n /3 operations 
where n is the number of unknowns. The solution of the Gaussian 
method involves eliminating the first unknown from the first 
equation and inserting it into all (n-1) remaining equations. The 
procedure will be described in an example of three equations with 
three unknowns using the modified Gaussian method. In the original 
Gaussian elimination, the Xj^ is eliminated from the first row, and 
according to the modified method term a^ x^^ is eliminated. So, we 
have 


20 



r 


a ll a 12 a 13 


X 1 


b l 

a 21 a 22 a 23 


X 2 

= 

b 2 

a 31 a 32 a 33 


X 3 


b 3 


which can be rewritten as 


a ll 

X 1 

+ 

a 12 

X 2 

+ 

a 13 

x 3 

rH 

A 

II 


a 21 

X 1 

+ 

a 22 

X 2 

+ 

a 23 

X 3 

b 2 


a 31 

X 1 

+ 

a 32 

X 2 

+ 

a 33 

X 3 

= b 3 

# 


Division of the terms of the first row with a^ yields: 


a ll 

a 12 

a 13 


a ll 

a l2 

a 13 

a 21 

a ll 

a 22 

a 23 


a 21 

a 22 

a 23 

a 31 

a ll 

a 32 

a 33 


1 

a 31 

a 32 

a 33 


(17) 


Elimination of a^x^ from the first equation gives 


a ll X 1 b l " a 12 X 2 ” a l3 x 3 


By insertion that into the second and the third equation we have; 


a lx * b l” a 12 x 2 “ a l 3 x 3^ + a 22 x 2 + a 23 X 3 b 2 


iff ( V a 12 x 2 - a 13 x 3> + a 32 x 2 + a 33 x 3 = b 3 ' 
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and, after reduction. 


a 21 a 2l a 21 

(a 22 - ^77 a 12 >x 2 + (a 23 - 377 a 13> x 3 “ b 2 ~ b l 


(a — — a 31 a 31 

32 a ll a 12 )x 2 + (a 33 “ a^ a 13^ X 3 = b 3 " a^ b l 


or, in a shorter form 


( 18 ) 


a 22 X 2 + a 23 x 3 


= b: 


1 1 , 1 

a 32 X 2 a 33 x 3 a 3 


The described procedure eliminates the unknown x^, i.e. the 
order of the system is reduced from n to (n-1) . Continuing 
this procedure by eliminating X 2 gives: 


(a 


33 


1 

1 

a 32 1 , ,1 

— a 23 )x 3 = b 3 

a 32 ,1 
1 b 2 ' 

a 22 

a 22 


or 


a 33 x 3 


= b. 


After the elimination has ended, the result is: 


i-W 
1— 1 

(d 

x i + 

a 12 

x 2 + 

a l 3 

x 3 = bj 

0 

+ 

1 

a 22 

x 2 + 

1 

a 23 

x 3 - b 2 

0 

+ 

0 

+ 

2 

a 33 

x 3 b 3 ‘ 


(19. a) 
(19. b) 
(19. c) 


The unknowns may be easily calculated from term (19) . From 
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equation (19.c): 


2 , 2 

x 3 to 3 /a 33 * 

Inserting x^ into equation (19.b) gives : 

X 2 = ^ b 2 ” a 23 x 3^ a 22 • 


Next; X 2 and x^ are inserted into equation (19.a) and x^ is 
calculated. 


Note that the elimination is done in the following way : 





which gives: 


a ll a 12 a l 3 


X 1 


b l 

0 a 22 a 23 


X 2 

= 

b 2 

0 a 32 a 33 


X 3 


b 3 


( 20 ) 


The elimination is repeated until the last unknown is reached. 

The important information is the number of operations required 
for the Gaussian elimination. 


According to the described procedure the following is needed 
( D ; division, I-I : multiplication, and S : subtraction ) : 




ELIMINATION 

FIRST 

REDUCED 

VECTOR 


SOLUTION 



COLUMN 

MATRIX 


B 


X 




D 

M 

S 

M 

S 

M 

S 

D 

X 1 

(n-1) 

(n-2) 2 

(n-1) 2 

(n-1) 

(n-1) 

(n-1) 

(n-1) 

1 

x 2 

(n-2) 

• 

(n-2) 2 

• 

(n-2) 2 

• 

(n-2) 2 

• 

(n-1) 

« 

(n-2) 

• 

(n-2) 

1 

• 

L 

• 

• 

1 

• 

• 

1 

• 

• 

1 

• 

• 

1 

• 

• 

1 

• 

• 

1 

• 

• 

1 

• 

• 

1 

x n 

— 

— 

— 

— 

— 

— 

- 

1 


n-1 

n-1 

n-l 0 

n-1 

n-1 

n-1 

n-1 



£i 2 

£i“ 

£i 2 

£i 

£i 

£i 

£i 

n 


i=l 

i=l 

i=l 

i=l 

i=l 

i=l 

i=l 


APPROXIMATELY 

n 2 /2 

n 2 /3 

3,_ 
n /3 

n 2 /2 

n 2 /2 

n 2 /2 

n 2 /2 

0 


The number of operations is defined with two types of sums of 

2 

real numbers , £ i and £ i j 


ii ) xi 

£i = - — — because the given sum is the re^l number 

arithmetic series which gives ^ , and 

is approximated with . 



i=l 


is the sum of the real numbers' squares, and is 
equal to 



\ 
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n 


1 ) 


which is approximated with y- . 
Thus reduction of matrix A requires 


n „ . n 
3 2 


D + 


n 


where only the first term y- M is kept, because it needs the 
most operations and the longest execution time. Vector b 
rearrangement requires 


T t! + I* s 

2 

operations of which only the first term, M, is kept. To 

2 ^ 

obtain the solution approximately y- M operations are needed. 
For example, a system of 30 unknowns will need : 

3 - 

y- + n = 9000 + 900 = 9900 

3 

operations. It is obvious that y- multiplications require most 
of the time and effort in calculations. 

A numerical example for the Gaussian elimination is 


2 2 3 


X 1 


1 

4 5 8 


X 2 

II 

2 

6 2 2 


X 3 


4 


1) 


n 


. 2 


An approximate term for E i 

i=l 


can be obtained from the integral 


% 

1 


2 di = d 

3 


n_ 

3 


_1 

3 


25 



2 2 3 


_ x i_ 


1 

2 ] 1 2 


X 2 

li 

i 

I o 
1 

l 

i 

m 


X 3 


X 1 


2 2 3 


X 1 


1 

— s 





2 N V 1 2 
X. 


X 2 


0 

\ 





3 -4 ; 1 


X 3 


1 


The solution is : 

x 3 = 1, 

x 0 = 0 - 2*1 - -2 f 

x x = (1-3.1 - 2 (-2) )/2 = 1. 


Checking is done bv inserting the solutions into the given 
equations s 

2-4 +3 = 1 

4 - 10 + S = 2 
6-4 +2 = 4. 

Regarding term (16), x = A * b, if vector b changes, the matrix 

A remains unaltered and the solution x is obtained bv a simple 

** -1 

procedure of multiplying matrix A, with the nev; vector b. When 

applying Gaussian elimination, instead of the inverse matrix A, 1 

the unknown vector x is calculated from term (19) . Though it 

appears at first that the vector b change requires the repeating 

of the whole Gaussian elimination, i.e. approximately n /3 

operations, this is not necessary ; and the number of operations 
2 

is only n (if only vector b is changed) . It is obvious that the 
coefficients a. . do not have to be calculated again in term (19) 
but only vector b (b. , b*, b^, ... b n )$ and for vector b 
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processing the terms in the columns which are placed under the 
diagonal are used. These terras disappear during tire elimination 
and zeros appear in their places. 

However, if we keep these terms, the calculation, because of the 

n 2 2 
vector b change, requires ■=— for solving vector b, and n /2 

^ 2 

operations for getting the solution, totally n operations. 

This Gaussian elimination supplement, by which all terms under 
the main diagonal are kept, leads te the so called LU 
decomposition. Let the matrix of term (19) be the upper 
triangular matrix : 


a ll 

a l2 

a l 3 


U 11 

U 12 

U 1 3 

0 

1 

a 22 

1 

a 23 

= 

0 

U 22 

u 23 

0 

0 

2 

a 33 


0 

0 

U 33 


The terms which appear under the diagonal during the elimination 
are called the lower triangular matrix L s 


1 

0 

0 


X 11 

0 

0 

1 

a 21 

1 

0 

= 

1 21 

1 22 

0 

•Ji 

a i 2 

1 


1 31 

1 32 

1 33 


where diagonal terms which are all equal 1 are added to the terms 
under the diagonal. 

It can be shov/n that : 


L U 

/vx 


A 




Substituting matrix with L U we proceed: 

A x = b , 


(21. a) 
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LUx 

Cs/ /V 


(21. b) 


which gives 


“ b , 



v — L ^ b • 

«L- r-f — 


Now , by inserting ^ into (21.6) we obtain 


x = 



Z • 


(21. c) 


(21. d) 


(21. e) 


Further examination of terms (21) reveals the following : 
term (21. a), i.e. matrices L and are obtained by the normal 
Gaussian elimination procedure? term (21. d) , i.e. vector is 
the symbolic representation of the vector b elimination procedure, 
term (21. e) is a symbolic representation of solution x from the 
reduced matrix A. 


3.2. Improved Precision of the Solution 

To increase accuracy, the obtained system A x = b solution is 
called x^ + x n » Inserting x^ + x n into the original equations 
gives: 


11 

*1 

+ 

a l2 

*2 

+ • • • • + 

a ln 


- 5 i 


b l 

21 

*1 

+ 

a 22 

*2 

+ .... + 

a 2n 

X 

n 

= b 2 

* 

b 2 

• 






. 


. 


. 

• 






• 


. 


. 

• 






• 


. 


. 

hi 

*1 

+ 

a ^ 
n2 

*2 

+ • • • • + 

a nn 

X 

n 

= b 

n 

* 

b 

n 


By this procedure the calculation accuracy is estimated, if all 

x. to x are correctly calculated, then all b's must equal b. 

in _ _ 

From the difference (b^ - b^) to (b n - b n ) , the calculation 
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precision is estimated. This simple control procedure allows 
improvement of the solution accuracy by subtracting the calcula- 
ted equation from the original ones 


a ll^ x l “ x l* + a l2* x 2 ” x 2^ + "*• + a ln* x n “ x n* ~ b l “ ^1 


a nl (x l " x l> + a n2 <x 2 ' x 2> + + a „n (x n " x n> “ b n " b n 


or 


a ll Ax l + a 12 Ax 2 + .... + a ln Ax n = Abj_ 


a , Ax. + a ~ Ax~ + .... + a„„ Ax, = Ab . 
nl l nz /. nn n ri 


Since matrix A has remained the same, and only term b, t b has 

1 n 2 

changed (into b.^ - b n ) the calculation is repeated with only n 
operations to solve for the unknowns Ax^ t Ax n » New values of x 
are: 


x l = x l + Ax^ 


x = 
n 


x + Ax 
n 


n 


This procedure can be repeated several times to obtain the best 
accuracy. 


3.3. The Inverse Matrix 

Using Cramer's Rule to obtain the inverse matrix is a very 

complicated procedure. The Gaussian elimination method requires 
3 

approximately n operations and is much more efficient. 
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It is done as follows : 

the system A x = b is given; the symbolic solution is x = a”^ b. 

The inverse matrix A, * elements will be called aij (to differ 

from the matrix A elements a. .). 

^ ID 

Supposing that matrix A” 1 is known, then : 


X 1 = a ll 

b l 

+ 

a l2 

b 2 

"t" • • • • 

+ 

a ln b n 

X 2 = a 2l 

b l 

+ 

a 22 

b 2 

• • • • 

+ 

a 2n b n 


( 22 ) 


n 


= a , b. + 
nl 1 


a n2 b 2 


a. b 

nn n 


If all vector b elements are zeros, except the first one which 
is b^ = 1, then the upper system is reduced to: 


X 1 a ll 


X 2 a 21 


(23) 


x = a , 
n nl 

If all b. =0 except b 0 = 1, then : 

l 2 

X 1 = a 21 


x 


2 


= a 


22 


X n = a n2 

Now, Gaussian elimination v/ill be applied to the above system 
A x = b, where b is : 
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The first solution values, x^ f x , have to be indentical with 
equations (23) , the second series values have to be identical 
with equations (24) , etc. The obtained solutions (x^ ~ x n ) there- 
fore represent the inverse matrix coefficients. The number of 
needed operations is : 

3 

- matrix A solution requires n /3 operations 

2 

- solving the right-side column, b, n /2 operations 

2 

- calculating x^ x r , n /2 operations 

Nov; the matrix A solution is not necessary any more. Continuing 
from the second column : 

- solving the second right side column b, but starting 

2 *” 

from bj = 1, because b^ = 0, (n-1) /2 operations 

- calculating x^ ~ x n , n/2 operations 

- solving the third right side column b, starting 

2 

from b- = 1, (n-2) /2 operations 

J 2 

- calculating x^ t x , n /2 operations, etc. 
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Totally, 


we need 


3 

n 

3 


+ 



(n-1) 2 . (n-2) 2 

2 + 2 + 



+ n 


2 

n 


2 


The above middle term is equal to : 


1 _ 

2 


n 

E 

i=l 


. 2 
i 




The total number of operations is therefore: 


n 


3 


T - 


+ 



n 


This requires only three times more calculations than is needed 
for the system A x = b solution. 


3.4. Gaussian Elinination Example 

To illustrate the simplicity of the Gaussian elimination the 
following FORTRAN program for solving the system of unknowns n 
is given; 
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n n n non non non non noon 


nl FOR. IS 

RESERVING FIELD FOR THE MATRIX A COEFFICIENTS. 
VECTOR 8 AND SOLUTIONS X 

DIMENSION A ( 100 * 100 ) » B( 100) .X( 100) 

READING MATRIX A C0EFICIENTS AND VECTOR B 

READ ( 5 * 100 ) N 

READ (5. 101) ( < A( I » J) ♦ I»1,N) »J*1*N) . <B< I ) *I*1.N) 

100 FORMAT (15) 

101 FORMAT (10F8.1) 

DO 15 J=1»N 
K=J+1 

SETTLING MATRIX A FIRST COLUMN AND VECTOR B 

DO 16 I*K.N 
A ( I » «J ) = A ( I . J ) / A ( J, J ) 

16 B( I )«B( I )-A< I.J)*B( J) 

SETTLING REDUCED (N-1)*(N-1) MATRIX A PART 

DO 17 L*K.N 
DO 18 M*K.N 

18 A(L.M)*A(L«M)-A(L. J ) * C J.M) 

17 CONTINUE 
15 CONTINUE 

CALCULATING UNKNOWN VECTOR X 
C*0. 

DO 19 J-l.N 
L=N— J+l 
K=N-J 

X(L)*(B(L)-C)/A(L*L) 

C-0. 

IF(K-l) 1.2.2 
2 DO 20 M*L.N 
20 C*C+A (K*M)*X<M) 

19 CONTINUE 

WRITING SOLUTION. VECTOR X VALUES 

1 WR I TE ( 6 * 102 ) (X(J).J-l.N) 

102 FORMAT (1X.10F12.3./) 

STOP 
END 
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3.5. Sparse Matrices 


Matrix G, which describes an electrical network has a great 
percentage of matrix coefficients equal to zero. This large 
number of zeros should be eliminated so that the operation is 
done only with coefficients other than zero. The small number 
of non-zeros (NZ) is the result of the electrical network na- 
ture. Matrix G is constructed as follows : 

4 



In one row there is a diagonal term and as many terms as the 
number of other nodes with which the observed node is connected 
(we do not count a connection to the referent node) . 

In the given example, the fifth row of matrix £ has 1+3=4 
elements. Electrical networks generally have 2 to 6 branches 
connected to one node. Suppose that the average is 5 elements 
per row and 50 nodes. Then, matrix G, will have 50 x 5 = 250 NZ 
elements and 2250 Z(zero) elements^ i.e.of 2500 elements there 
will be 10% NZ and 90% Z. It is clear that such a great number 
of Z elements should be eliminated to obtain a more efficient 
calculation. 

The second property is that matrix £ is symmetrical regarding 
the main diagonal. Since matrix £ does not contain any dependent 
current source, it is symmetrical. However, if the network 
contains dependent current sources, then there are some non- 
symmetrical terms too. This may be demonstrated by the Ebers- 



Moll's transistor model (Fig. 10) . 



E 11 


Fig. 10. Ebers-Moll 3 s transistor model gives a 
symmetrical structure to the G matrix. 


In the given example every dependent source gives three non- 
symmetrical terms. It is fortunate that branch 5,7 and branch 
5,11 are interdependent. In such cases, where there is a recipro- 
cal dependence of two branches, the matrix is symmetrical; however, 
symmetrical coefficients do not have the same values. Thus when 
considering sparse matrices their symmetrical structure will be 
taken into account. 


Thus there are over 90% zero elements, matrix C has symmetrical 
structure, and the number of non zeros is different per given row . 
Gaussian elimination should be performed in a way such that it 
does not occupy too much of the memory, and that it does not 
require too many operations. Furthermore note the following s the 
matrix Q coefficient g^^ will be processed only when frontal coef- 
ficients in the row i, and column g , differ from zero,, regardless 


of whether the coefficient g. . 4 0, or 

-’l-j r r 

after processing, the element becomes g 


g • -=0 . 
ID 

ij * 0 


In case g . .=0 
D 

so a new non- 


zero element is generated. The fact that Gaussian elimination 


35 



generates new non- zeros leads to the growing number of elements 
in matrix <3 during the elimination. The purpose is therefore to 
arrange rows and columns so that a minimum of non-zeros is 
generated during the elimination. This will occupy the least 
amount of words in memory, and will require the least number of 
operations. 

Since the Gaussian elimination procedure has been covered, stor- 
ing and arranging of the matrix G elements will be demonstrated. 

a) Matrix G element storing 
Non-zero elements are read by rows as follows : 

A ( 1 ) = x 
A (2) = x 


A(8) = x 


-zeros, then 

n words are needed for the row index, R 
NZ words are needed for the column index, C 
NZ words are needed for the element value, A 

i.e. totaly n + 2NZ words. According to the given example with 


row 


column value 


1 

1 

X 


3 

X 


7 

X 

2 

2 

X 


3 

X 


7 

X 


9 

X 

3 

5 

X 

• 

• 

• 

• 

• 

• 

• 




R (1 ) = 1 C(l) = 



If matrix G is n x n, and NZ is the number of no 
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n = 50 and 10% NZ elements, this gives 50 + 250 + 250 = 550 
words (instead of 25oo) 

b) Matrix G arrangement 

In general the matrix G is arranged so that the row and column 
with the least elements is in the first row. The row and column 
with the most elements go in the last row. Rows and columns are 
used because in the matrix with symmetrical structure, row g and 
column g have the same number of elements. 

Therefore, switching the row means that the column must also be 
switched to keep matrix G structurally symmetrical. 

After this first arrangement is done, there are several methods 
on how to proceed with the G matrix rearrangement, one of which 
is given in reference (5) . Let us see one of the possible meth- 
ods whose principle is given in reference (6) . 

The top of the matrix G, has row with the least elements, say it 
has one off-diagonal element. This r ow does not generate new 
non-zeros. Regarding Fig. 11, after eliminating the first row 
and column, row k will lose 1 non-zero and will become a row 
with 1 off-diagonal element too. Thus row k , too, should be 
switched to the top of the matrix. In other words, after the 
first unknown is eliminated we treat the rest of the matrix G, 

.V 

i.e. (n-1) (n-1) matrix as follows: we scan (n-1) rows left 
and the row with the least elements is put on the top, and 
the second unknown is eliminated. If some new elements according 
to the elimination nature are generated they should be added 
into the proper rows. The same procedure is repeated for 
(n-2) (n-2) matrix and so on. At the end a row elimination 
order is obtained which has to be stored in a separate file. 

This file is used when the Gaussian elimination of the same 
matrix £ is to be repeated. This procedure gives very good resu- 
lts with respect to minimum number of operations and minimum number 
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of new generated elements (so called "fill-in”) . The very minimum 
is not obtained and as far as is known such method has not yet 
been proposed at present. 


“lk 


x 


X 


X 


1 


kl 



X X 

X X 

X 

X 


Fig. 11. An example showing that vow and column 
with two elements do not generate a new 

element . 


An example for matrix with 24 unknowns and an average of 4 elem- 
ents per row (including diagonal term) is given. Matrix is gene- 
rated randomly using the Monte Carlo Method (Fig. 12. a). Gaussi- 
an elimination result executed on the originally given matrix 
is shown on Fig. 12. b and the number of multiplications (and 
subtractions) is 452, number of division is 86 and number of 
fill-ins is 100. VThen the reordering scheme described above was 
applied, Gaussian elimination result is shown in Fig 12. c and the 
number of multiplications (and subtractions) is 131, number of 
divisions is 51 and number of fill-ins is 30. On the left column 
of Fig. 12.c the row sequence elimination is given. 
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-8e 


123*5678901 2345 6789 01 234 
IX. X 

2. X 

3. .X.X X-.. 

*» •• X 

5. »X • X X • • • .X..*X ..XX ...X* 

6. •• .XX....X..*. ...X 

7 X 

8 X....X 

9 XX.. .X.X 

10 .xx x ... . 

11. X..XX....X....X 

1 2. «X ........ Xw*. *X ...... « 

13. X....X-X 

14X x • . • • XX • .-x • . «x • 

15. XXX 

16. XXX....X 

17 .XX. XXX.. 

18.. ..X.X. XXX- *X.. 

19. . ..XX.X.....X* ..XX.X*.. 


20. X..X-... 

21. .XX.«... X-X.X.*. 

22 XX...X.. 

23....X........X XX 

24 XX 


a) 


123455789012345678801234 


lx 

2.X X 

3. . / .y ...... x. x... 

4. .. X X... 

5.. X.XX*.* . XX. -X.. XX. X«X. 

6. . . .XX* . . .XX. .x .*xx • x • x • 

7 

8 • «X X 

9. . ..... .XX. ».X. X. ...... . 

10 xx...x.y...x.... 

11. x» .XX... .XX. .XX. XV. X.X. 

12. . x *xx. •• *xx. .xxxxx «x *x« 

13 x. . . .x .x .. .x 

1 4x » » .» .. »yx. ..xxx»* xx». x * 

15. . . .xx... .xxxxxxxxxxx.x. 

15 XXXX.XXXXXXXX.X. 

1 7. . ....... ..x. .XXX XX xx xx. 

1 9. .. .XXX.. .XX. -xxxxxxxxx. 

1 9 . . . .XX.X. .xxxxx XX XX XX xx. 

20 x. ..xxxxxxxxxx. 

21. .XXX x. .. .xx. .XXX XX XX xx. 
22 XXXXX xx. 

2 3. .. .XX.... XX. XX XX XX XX XXX 

24 XX 


b) 


124 74 82 30 233654 681799015 


IX 

2.X 

*..X. 

7...X X 

24.. • 

8 X.X... «... 

12 X....X X..... 

13. . ...X.X X*.. 

20 .. . 

22 X.X 

23.. ..X.....X...X. .X 

3 X....X XX. ...X 

S. x-xx 

15. X.....XX X..X 

14X.........X*. XX . • . • XX • • X 

16. 

18.. X X.X-X..X 

21.. X........X*.. .* XX • X • • X 

17. XX.X.... XXX- XX.X 

XX- ..XX XXX 

is. •••* x»x. ... XXX .xxxxx XXX 

10 X......X..XXXXXX 

l i« x> ......... x« * x... XX XX X 

s -XXXXX XXX 


a) 


Fig . 12. a) Randomly generated 24x24 matrix with an average of 4 
nonzero elements per row y 

b) Gaussian elimination result executed on the given matrix 3 

c) Gaussian elimination result after reordering . 
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4. SENSITIVITY COEFFICIENTS 


After the node voltages are defined, it is often necessary to 
calculate the sensitivity of an individual voltage to the network 
parameter changes. A sensitivity analysis is usually done for on- 
ly one node, i.e. the output voltage. Here, the described proce- 
dure concerns the linear networks and is designed for use of 
the total differential. Let V o = f(x^, Y^r ••• x n ) be an output 
voltage; x^ supplies the network parameters - resistors and 
current sources. The total differential from V q is 


dV o = M7 dx l + Hq dx 2 + 



The solution requires finding the differential dV Q , but from the 
matrix term for the network. 


The network is described by the term 


A A V = A I 
^ 'n/D ^ — n — s 


and the total differential is equal to 


A dG„ A fc V + A G„ A t dV = A dl 

^ ' — n ~ D — n — ' — s 

which, with G = A G_ A fc , follows for dV 

dV = G -1 A dl - G -1 A dG^ A fc V , 
— n ^ ^ — s ^ -'D ~ — n 


It 

G -1 A 


dl - dG„ A t V„ 

1 



— s D — n 


Let the matrix G * A = H and V = V, . The upper term is written 

r*s — n — U 

in the extended form (n is the number of nodes, and b is the 
number of branches) as: 
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dV. 


1 


• 


• 


dV 

■ = 

r 


dV 


n 



h ll * * * h lj •••* h lb 


h rl * * * h rl •••• h rb 


h nl * * * h nj ***• h nb 


dI sl - V bl dG l 
dI 8j ■ V bj dG j 


dI sb - V bb dG b 


Matrix H represents the sensitivity coefficients matrix, because 


dV r 

’i . dl . h rj * 

rj S3 J 


S G . " HG 1 " _h rj V bj * 
rj 3 j j 

S is the voltage sensitivity of the the r-th node in relation 

rj 

to a current source in the branch j. S is the voltage sensiti- 

G rj 

vity of the r-th node in the relation to the branch j conductivity. 

To find the matrix H = G *A, since matrices G and A have been done, 

the only problem is the matrix J3 inversion which requires 

3 /s/ -1 

approximately n operations. The multiplying of matrix G by 

A and other operations is not a big burden for the computation. 

To determine the sensitivity coefficients of one output voltage 

(for example node voltage) , there is no need to calculate 

the entire coefficient matrix H', but only its row r, because: 


dV = h , 
r rl 


< dI sl- V bl 


dG 1 )+h r2 


(dI s2- V b2 dG 2» + - 


• +h rb (dI sb- V bb 


acy 


Before the above mentioned H matrix r-th row is obtained, the 

_i 

r-th rov; of the matrix G must be done. The previously described 
matrix inversion procedure, using Gaussian elimination (see page 
29) , gives one column per one inversion step. To calculate the 
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sensitivity coefficients of one row, the matrix G inversion should 
be executed by rows* This is possible and is done with the help 
of the following relation: 


(G t )" 1 = (G _1 ) t 


The purpose of this relation is to invert the transposed matrix 
G \ by the usual procedure, i.e. columns. After this, the inversed 
matrix G* 1 column is equal to the inversed matrix G row. Thus the 

r* v 

sensitivity coefficients of the desired output node voltage are 
n 3 

obtained by -j— operations. 

The sensitivity application and its connection with the relative 
Orror of a very simple example is as follows (Fig. 13) : 

v s = 10 V, 

D 

1 R = 10 Q G 2 = 0.1 l/Q 

— v i R 2 = 5 0 , G 2 = 0.2 1/fi 

R 2 

resistor tolerance - 10%. 

Fig. 13 . Sensitivity computation example. 



Voltage sensitivity to the parameter is 


S 


G 


11 



t 


A - Gi +G - 


1 =io SLil = 3,33 V 


0.3 


dV, 


G. 


G n dG i ° (g 1 +g 2 ) 2 


= 10 


o.2 


(0.3) 


= 22,2 . 


If the resistor R^ changes for + 10%, then G^ changes from 0.1 
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to 0.091 so G^ = 0.009 1/J2 , and the output voltage difference 
is 

AV, = AG, = 22.2 0.009 = 0.202 V 

1 G 11 1 

Working only with the relative errors (tolerances) , the term is 
arranged like this : 


iV l . S G , G G 1 . AG 1 
V 1 " V 1 1 G 1 " G 1 


, fl 
’ G 1 V 1 


where 


AV, 


AG, 


v P VI' and Gl P G1 


VI 

which gives 


Gl V, 


Gl ' 


VI 


22,2 


0.1 

3,33 


P G1 " °‘ 667 P G1 


This means that if R changes by 10%, the voltage V^ will change 
by 6.06%, and in total this amount gives V^P V ^=3 . 33* 0 . 0606=0 . 202 V. 

The observed change of value V^ and the relative change P v ^ 
are obtained from the same sensitivity coefficient, the only 
difference is in the form in which the solution is given. 
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5. NONLINEAR DC NETWORKS 


The nonlinear network, other than the linear elements, contains 
nonlinear elements such as diodes and transistors. A nonlinear 
element cannot be directly presented in matrix g with a constant 
term. So, in the nonlinear network the operating point of some 
nonlinear element is defined by Newton-Raphson's method. The meth- 
od consists of the following: to find an intersection of func- 
tion y = f (x) (Fig. 14) with the abcissa, an x^ is chosen and 
y ^ = f(x^) is defined (point A). 



Fig. 14. Example of Newton-Raphson* s 

method. 

A tangent y'=f'(x 1 ), is drawn through point A. The intersection 
of the tangent and abcissa is X 2 » The procedure is repeated un- 
til either y or the absolute difference (X 2 - x^) = Ax becomes 
as small as is required. 

The procedure is now adjusted for analyzing nonlinear electrical 
networks with a larger number of elements. It is understandable 
that in a complicated network, gradually approaching the solution 
cannot be illustrated graphically. The entire procedure will be 
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With a known diode voltage. 


V JD 


kT 


I D = I s (e V T - 1) , V T = = 26 mV 


dl r 


v 


G n = WT = (I - /V -> e V m = dr. + IeJ/Vn, . 


S' T 


T 


Through the point Ip, Vp a line with the slope G p is drawn: 


(I " I D ) " G D (V “ V ' 

I — G (V - V) + I„, if V=0 we have 


*00 - X D - G D V D 


D' 

or 


Z D 0 - S * (V T - V - *s - g d v t (1 - ln TT’^S 


Between nodes m and n, instead of the diode, the conductance G D 
and the current source Ip Q are connected. 

In the example of Fig. -16, the diode initial values are I D , V D , 
(Fig. 17. a) . 


The line which corresponds to conductance Gp intersects the 
conductance line, G, at voltage Vp. The current is calculated 
by Ip = f (Vp) , a tangent is drawn through the new point on the 
diode's characteristic, and the intersection with line G is 
found. 

The procedure is repeated until voltage AV = (Vp - V Q ) is within 
the selected error limits. Voltage AV is the tolerance (usually 
50 ^iV) because it is not possible to calculate the voltage Vp 
exactly. 
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a) b) 

Fig. 17. Newton-Raphson* s method application 

example . 


In a complex network, in place of each nonlinear element, the 
conductance G D is placed parallel to the current source I DQ . The 
system is solved using the node voltages t V . From the node 
voltages the nonlinear element's (branch) voltages are calculated 
using V . After the value V^, i.e. V D , for every single 

nonlinear element is obtained, tolerance voltage AV=(V^ - V^) is 
calculated and, if it is within the 50 pV tolerance the procedu- 
re is completed. If not, voltage V D for each nonlinear element 
whose difference voltage AV is not whit in the given tollerance is 
set equal to voltage V£ and new values G Q and I^q are calculated 
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and inserted into the network. The procedure is repeated until 
the difference voltage Av for every single nonlinear element is 
less than 50 fCJ . 


Electrical networks nonlinear elements are diodes and transistors 
that are modeled using diodes. Because of this we have to notify 
some exponential function e V D/v t characteristics. During the 
elimination, some difficulties can arise in the calculations if 
and V D (new and old diode voltage values) differ by more than 
2 V T , i.e. (V^ - V D ) > 2V T .This is explained in the following 
example : 


*D " Z S e V T ' 


% v d- | ' 2v t 

X D = I g e V T X s e V T 


Then I d/Id = = 7 >8 • 

Due to the exponential dependency, a new current, 1^, is eight 
times larger than the current 1^. Such a large difference is un- 
desirable and therefore the growth of the diode voltage is 
restricted to 2VT = 52 mV, i.e. = V Q + 2V T (instead of the 
calculated voltage V^) 

A modification of Newton-Raphson's method provides better answers 
to the pecularities of the diode characteristics. This is shown 
in Fig. 17. b. 

After obtaining the solution the nonlinear element current is 
calculated using 

— ~ <£d — b - ~D — n' i,e * ip and, using these currents, 

the new values G D and I DQ are calculated. By this procedure, the 
solution with fewer iterations is obtained. 
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6. LINEAR AC NETWORKS 

Before beginning linear AC network analysis explanation of how 
any nonlinear network may be linearized must be given. Generally, 
networks are nonlinear because, in addition to R, L, C, elements 
they also contain nonlinear and active elements ) diodes and tran- 
sistors. In this type of analysis, nonlinear elements must be 
linearized about the operating point and the values inserted 
into the network. In the diode linearization example, the tangent 
of the operating point according to the known term is calculated 
as shown : 


G 


D 




V. 
e Vi 


D 


(I D + V /V T ' 



Fig. 18. Diode linearization for AC 
analysis. 


DC analysis of a nonlinear network requires the diode model 
according to Fig. 16, while AC analysis needs only the conductance, 
G d (Fig. 18) . Once we have performed DC analysis, all models of 
nonlinear elements needed for AC analysis are existing. 
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After the DC state is defined, all nonlinear elements are 
introduced into matrix G with adecuate conductances. G_, and 
current sources, I DO * Conductances, G^ , are kept in the network 
while current sources I DQ are removed. After this procedure, 
linear AC analysis may begin. 

Generally the matrix terms of DC analysis are used. This provides 
the admittance matrix Y , instead of the conductance matrix Cj, and 
voltage V n and current sources I_ vectors are complex. So, we 
write : 

Y V* = I* - 1} (27) 

— n — sn 


where 


Y = 


G. . + j B. . 
3-D 3-D 


v* = V„ + j V 
— n — n J — n 


I* = I + j I 
— sn — sn J — sn 


The upper term could be used as it is, in which case the first 
row of term (27) will be : 


< G ii + i B n> ( v i + 3 V + (g 12 +j b 12 hv 2+ 3 v 2 ) + 


* * ^snl + -^ ^snl 


( 28 ) 


This form is inconvenient for the work on the computer, so it 
must be rearrenged. First, the standard branch for AC networks 
is defined (Fig. 19) 


Mark V* represents the voltage having the real component V and 
the imaginary component V. 
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I 


+ 

s 




Fig. 19 . A standard branch for AC 
networks . 

A diagonal admittance matrix is 


Y 

-vy 


D 


G. . + j B. . 
11 J 11 


and can be presented as the sum of matrices G D and B D as 
follows : 

Y = G_ + j B n • 

Now, the procedure for obtaining the admittance matrix can be 
applied usnig the diagonal admittance matrix, i.e. 


Y = A Y n A t = A (G + j B) = 

= A G^ A fc + jAB^ A t = G + j B . 
^ < D - — ^ 4 — ' 


(29) 


This term implies that the admittance matrix can be constructed 
from two matrices, and B, where each one is calculated independen- 
tly using the incidence matrix A. Thus, matrix G, is constructed 
using the previously established rule (see page 8) without regard 
to reactive branches of the network 0 Matrix B, is constructed using 
the same rule without regard to resistive branches of the 
network. Using this procedure, term (27) can be written as: 
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(30) 


w + 3 B) (V n + j V n ) = <I sn + j I sn ) , 


or 


G V 


V + j B V + j G V - BV = I 
— n J — n — n - — — n — : 


+ j I 

sn J — sn 


(31) 


The real part is equalized with the real one and the imaginary 
with the imaginary one: 

G V - B V = I 

— n — n — sn 

B V + G V = I 

/N/ — n — n —sn 


Which gives the matrix term: 



(32) 


The upper term is used directly for solving the linear AC 
networks. The first row, written more extensively, reads: 


G 11 V l +G 12 V 2 


+... + G 


In 


V n - 


B 


11 


V 1 - 


B 12 V 2 


“ B ln V n 


= I 


snl 


(32. a) 

This form is much more suitable than that of term (28) . 

In term (32) note that the number of unknowns is two times 
larger than the number of nodes because of the fact that all node 
voltages have real and imaginary components. Matrix is four times 
larger than matrix because instead of (n x n) elements it con- 
tains (2n x 2n) elements. Thus, the calculation is also larger. 

Applying the described procedure on the network in Fig. 20. 
matrices G and B are written as: 
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Fig. 20. The linear AC network example. 
where B 2 = u C 2 , B 3 = gU , 

Ij^ - I 1Q sin tot , I 3 = I 3Q sin (tot + 4 >) . 
The entire matrix term which describes the network is: 
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or, for the first two rows: 


G 1 V 1 ’ G 1 V 2 


= I 


si 


= G 1 V 1 +(G 1 +G 2 )V 2 - (B 9 +B^) v 9 + B, V, - , etc. 


2 3 


2 3 s3 


Here, the analysis is performed at the chosen frequency. 

If we change the frequency, the matrix G, remains the same, while 
the matrix B changes. The terms of matrix B are of the form 
toC, ^ or their sum and they change with the frequency. 

This means that for every frequency value, matrix ^ elements 

must be recalculated and inserted into term (32) , and new values 

for node voltages V = V + j V are calculated. Regarding the 

fact that the observed frequency area is usually large, the 

frequencies are shown in logarithmical scale. The frequencies change 

as f =10 f . where x is calculated from the number of 
n o' 

observed frequencies per decade. For example, if 4 frequencies 
per decade in the same logarithmical increments are desired, 
then 4x = 1, x = 0*25 so 


f 

f 

f 

f 

f 


0 

1 
2 

3 

4 


10° f _ 
o 


= 1.000 

f o 

10°.25 

f o 

= 1.780 

f o 

10°* 50 

f o 

= 3.165 

f o 

io°. 75 

f o 

= 5.630 

f o 

10 1 * 

f o 

=10.000 

f o 


Input and output voltages and currents. 


(V+jV) and (I+jl) , can in 
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some programs be changed to \|v 2 + V 2 /$ v and l/l 2 + I 2 /±i . 

Thus by changing the input signal frequency as a result of AC 
analysis, the amplitude and phase network characteristics in 
any node are obtained. 

For dependent current sources, and voltage sources, the principle 
described for the DC analysis is used, except that consideration 
should be given to their complex nature. 
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7. TRANSIENT ANALYSIS OF NONLINEAR NETWORKS 


In this chapter a description is given of the procedure for 
analyzing the nonlinear electrical networks which, in addition 
to the resistance R, contain capacitance C, and inductance L. 

All three kinds of elements, R, L, C can be nonlinear. The input 
signal is the time function. Thus the transient analysis of 
nonlinear networks will be discussed. 

Up to this point, a static analysis of nonlinear networks has 
been examined. This means that only the resistor, linear or 
nonlinear, was included. Next, the networks with constant 
capacitances C and inductances L were included, as well as 
linearized nonlinear elements and the input signal was sinusoidal. 

First, the transient analysis of linear networks is described, 
the solution stability relating to that is examined, and the 
networks which contain the nonlinear C and L elements are 
analyzed. Finally, the nonlinear resistances are added to the 
network . 


7.1. Transient Analysis of Linear Networks 

Suppose that Fig. 21 represents one part of a given network. 





which, when inserted into the equation of the current balance in 
node 2, yields: 
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<'V V 2 ,dt - VVV 



- c. 


dt 


- V tJ 


This simple example shows that the old system of algebraic 
equations, G V n = I , will change to a system of integral - 
differential equations. Using the computer the analytical 
solution is unobtainable? this means that the integral-differential 
equation must be changed into an algebraic one and numerical 
method should be applied. The mathematical approach that was 
used for solving the resistance networks is also employed for 
solving dynamic (linear and nonlinear) electrical networks. 

Solving the differential equations is called integration. There 
are many numerical integration methods. Three of the most common 
are described here. When choosing the integration method, 

the one which is simple enough for the application and which 
gives a stable solution should be chosen. The first of the 
three methods is Euler's method, which is often unstable but it 
is the least complex. This method is used to describe the principle 
of the differential equation's numerical integration. 

For an example of Euler's method, refer to Fig. 22. 




Fig. 22. Simple RC circuit. 
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The mesh equation is: 


U = I R + 
s 


H 


I dt , or 


dU C 

U s " RC + U C 


(34. a) 


Equation (34. a) is the differential equation which connects the 

variable U c and its derivative dU c /dt at any time, t. If the 

signal voltage U and the voltage on the capacitor U_ are known, 

s u 

the voltage derivative U c can be calculated from equation (34. a) 
as follows: 


fV ^ ^ 

dt RC RC 


(34. b) 


Now, the variable U c value in time t+At is calculated using 
Euler's method: 

dU p (t) 

U c (t+At) = U c (t) + At — ^ — . (35) 

The voltage U c curve is approximated by the tangent. The 
approximation becomes more accurate as the time interval At 
decreases. After the calculation, the new voltage value U c 
is inserted in term (34,b) and the new derivative value is 
calculated. This continues until the end of the given time 
interval is reached. The error which appears according to term 
(35) can be calculated by using Taylor's series: 


U c (t+At) = U c (t) + At (t) 



u; 


(t) + 


u. 


’(t)+.. 


(36) 
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The error is, as is shown in the comparison of terms (35) and 
(36) , equal to: 


..3 

At , At S , 

2T U C (t) + TT U C (t) + 


At 2 

Thus it is proportional to -y U^"(t). To avoid forming differential 
equations, and to keep the mathematical procedure that was used 
for resistive networks, equation (35) is written as: 


dU 

T Vl - °n + At -a? ' 


U . . = U + I 
n+1 n C n 


(37) 


An equivalent circuit for the capacitance C following from equation 
(37) , is shown in Fig. 23 which is a voltage source. 




m 

<? 


© ] 


f(V i L Lu n) 





Fig. 23. a) An equivalent circuit for the capacitance C t 
b) An equivalent circuit for the inductance L t 
according to Euler* s formula. 
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A calculation of the example in Fig. 22 (using the principle shown 
in Fig. 23) is given in Fig. 24, 


R=10fl 




n At 

U (V) 
n 

(A > 

< U n + l> (V) 

0 

0.0 

1.0 

1.0 

1 

1.0 

0.9 

1.9 

2 

1.9 

0.81 

2.71 

3 

2.71 

0.729 

3.50 

4 

3.50 

• 

• 


Fig, 24. RC network equivalent circuit with the ■ 
numerical results. 


Note that the network containing capacitance C is now included 
in the network that contains only voltage sources and resistors. 
Mathematically, this means that the differential equation is now 
an algebraic one. 

For the inductance L, the fundamental relation between the voltage 
and current, U_ = L dI T /dt. Euler's integration will give 

Jj Xj 

dl 

I n+1 = I n + At “d? ' respectively 
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(38) 


I 


n+1 



An equivalent circuit for inductance L is, using tern (38), a 
current source and is shown in Fig. 23. b) . Because of the 
capacitance C, the described procedure. is not convenient for the 
nodal analysis application since it includes the voltage sources 
in the network. 

The second integration method is a modified Euler's method, called 
Backward Euler (B. E.), which is more convenient and gives a stable 
solution of the differential equation. Equation (35) , using this 
method, is written as: 


U c (t+At) = U c (t) + At 


dU c (t+At) 

dt 


0 . . = U + At U . . 
n+1 n n+1 


and, after substituting U^ +1 “ I +1 /C, gives 


U n + 1 - U n + ^ I n+1 * 


respectively 


(39) 


This procedure is justified because the relation between voltages 

U . , and U can be a derivative in time t , just as it can be a 
n+1 n n 

derivative in time t +1 » From term (39): 


I 


n+1 




U 

n 


(40) 


which leads to the equivalent circuit as in Fig. 25. a. The current 
I +1 is equal to the sum of currents through the conductance G c 
and the current I c , which satisfies equation (40) • 


I 
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m 



4 


m 



a ) 


b ) 


Fig. 25. a) An equivalent circuit of the capacitance C t 
b) An equivalent circuit of the inductance L t 
according to the Backward Euler method. 


Similarly, with inductance L: 


*n+l " + 4t J n + 1 


i = i + u . , 
n+1 n L n+1 


(41) 


which gives an equivalent circuit as in Fig. 25. b. The Backward 
Euler method is much more suitable because both C and L are 
presented with the same kind of equivalent circuit and they 
satisfy demands of the network calculation by the nodal method, 
i.e. there are no voltage sources. Because of this, the Backward 
Euler method is used in many programs for dynamical electrical 
network calculations. The error generated by this method is the 
same as that of the original Euler's method, i.e. it is 
proportional with 


62 



/ 


respectively 



(t) 



(t). 


The third method is the trapezoidal integration method which can 
be explained in two ways: 

a) The voltage of the capacitance C in time t+At is calculated 
by 


dU (t) dU r (t+At) 

u c (t+At) = U c (t> + if (-S4- + — S-35 ), 

respectively U n+1 - 0 n + if (IT' + U' +1 ) , ,«. a) 

1 I , , 

or, by substituting U' = , U' +1 = — 

U n+1 '°n + H (I n + W • (42 - b) 


This procedure uses the derivative mean arithmetical value in 
time t and t+At. From term (42. b) 


I = — — U - (— U + I ) 
n+1 At u n+l v At U n n' 


(43) 


which leads to the equivalent circuit in Fig. 26. a. The current 
I is equal to the sum of the currents through the conductance 
G c and the current I c , which satisfies equation (43) . 

Similarly, by manipulating the inductance L: 



m 

0 


m 




n 



b ) 


Fig. 26. a) The capacitance C equivalent circuit t 
b) The inductance L equivalent circuit 
according to the trapezoidal forumula. 


n+1 


I n + *5 (I n + ^ + 1» ' 


or substituting 


„ _ U n+1 

n L ' J ‘n+1 L 


I = I + (U + U . _ ) . 

n+1 n 2L n n+1 


(44) 


After arranging the right side terms ; 


t = AA n + (AA it + I ) 
n+1 2L U n+1 + l 2L U n V ' 


(45) 


which leads to the the equivalent circuit in Pig. 26. b. The current 
In+1 * s e< 3 ua ^ to sum t ^ ie currents thorough the conductance 

G t and the current I T , which satisfies equation (45) . 
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The error which appears in the described integration procedure 
(formula 4 2. a) is calculated by comparison using Taylor's series 


°„+l - °n + “5 K + °n + l> ' 

U n+1 = °n + 4t U n + T <°n + l " U n> 


(46. a) 


Now, using Euler's method the term for the first derivative in 
time t+At is given as 


U ' - U' + At V'' 

n+1 n n 


t 


from which 


U 


n+1 


“ U n " At U i 


inserting this into the 
term (4 6. a) gives 

A . 2 

U _. A1 = U + it U' + 2%- U''. (46. b) 

n+1 n n 2 n 

This integration method is more accurate than Euler's methods, 
because the function is approximated with the first three terms 
of Taylor's series, i.e. it includes the square term, so that the 
error is proportional with At /31 . 

b) The second approach to the trapezoidal integration gives 
the same formula. 

The voltage of the capacitance C, at ime t+At, is equal to 
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t+At 

U(t+At) = U(t) + ~ | I dt = 

t 

= U(t) + ^ (I (t+At) + I(t)) 


or , simplified. 


U 


n+1 


= U n + M (I n + W' 


(46. c) 


This is the same as term (42. b). In Eq. (46. c) the voltage growth 
on the capacitor is calculated by integrating the capacitor 
current by the trapezoidal rule (Fig. 27) . 



Fig. 27. Trapezoidal integration. 
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7.2. Integration Method Precision and Stability 


Precision and stability of the three described methods will now 
be explained by a simple example of the RC circuit (Fig. 28.) 



E = U + V 
t = nAt 


Fig. 28. Simple RC circuit 


a) Euler's method (E) 


At any moment, t = nAt, 


E = U + V , or U = E - V 

n n n n 


(47. a) 


and at t+At = (n+l)At 


E - °n + l + v n+l 


' or D n« - E ' v n+l- <47 ’ b) 


and also 


I 


n 


V 

n 

R 


(47. c) 
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The voltage of the capacitor is 


U 


n+1 


= U 


. At T 
+ C" T l 


(47. d) 


Using the first three terms (47. a,b and c), the correspoding 

values for U , U and I are inserted into term (47. d) and 
n+ i n n 

give: 


E ' V n + 1 “ E * V n + BC V n ' 

V n+1 _ V n U ‘ W* ( 

If the voltage, V , is to be time dependent, then (t= RC) : 

V 1 - v o (l - - E o « - ' 


because at t = 0, the starting voltage V Q on the resistor equals 
the source voltage E q . Then 


V 


2 


V, (1 



E o (1 



f 


and, applying this procedure n times gives 


At. 


n 


V = E o (1 - . 

no t 


b) Backward Euler method (B.E.) 


( 49 ) 
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Equations (4 7. a, b and c) are the same, except that the voltage 
of the capacitor is: 


Vl - D n + B* r n + l ' 


which after the insertion gives. 


E — V . - = E 
n+1 


V n + If V n + 1 ' 


V 


n+1 




(50) 


Thus the voltage at time n^t following the steps used in Euler's 
method is equal to 


V 


n 


E 


o 


1 


(1 


+ 


At > n 
t ' 


(51) 


c) Trapezoidal method (TR.) 

The equations (47. a, b and c) are the same, except that the 
capacitor voltage is 


U 


n+1 


= U n + M (I n + I n+1 ) 


which, after the insertion, gives 


69 



E - V = E — V + — V + ~ V 

n+1 n 2RC n 2RC n+1 9 


V . . = V 
n+1 n 


1 - 


1 + 


At 

It 

Tt 

7F 


( 52 ) 


Thus the voltage at time nAt following the steps used in Euler's 
method is equal to 


V 

n 



(53) 


The exact analytical solution for the observed RC circuit is: 


V(t) 


E o e 


-t/T 


i.e. substituting t = nAt 


V 


n 


E o e -n4t/T 


(54) 


Using equations from (48) to (54) an estimate of the accuracy and 
stability of all three described integration methods may be made. 
First, the precision which is obtained in the case when the time 
increases by At is evaluated. The results forAt = 0*lx and 
At =t are given in the following table: 
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PRECISION 


(i - £|) 


0.90000 


0.0 


0.90909 


0.5 


TR. 

exact 

1 - At/2r 

-At/t 

1 + At/2r 

C 

0.90476 

0.90483 

0.333 

0.368 


The trapezoidal method best approximates the exact solution. The 
stability can be estimated by increasing the time step At and 
observing the result by odd and even n's. The three methods are 
compared in the following table: 

STABILITY 


stable for 


for very 
large At 


(1 ~ &) 


0 < At < t 


( -£> 


divergent 
and osc. 



B.E. 


1 

(1 

* , n 

+ ££) 
t ' 

0 

< At 


stable 


, . At / 
1 + T r 7 


0 < At < 2t 


( - l) n 


stable 
but. osc. 

























Euler's formula is stable when the time step is At<T . When At 
is very large, the solution oscillates and diverges more as n 
increases. This method is called an "explicit" integration method. 
It is generally stable only if the condition At<r is satisfied. 

If the electrical network has more time constants, then time step 
At for the stable solution must be smaller than the circuit's 
smallest time constant, At<T min . 

The Backward Euler method is stable for any time step, At, and 
the solution's precision decreases as the coefficient At/t 
increases. This method is an "implicit" integration method which 
is generally stable. 

The trapezoidal formula is the most precise and is stable even if 
the numerical solution oscillates for At>>T. This method is also 
an "implicit" integration method. The trapezoidal formula is 
preferrable to the Backward Euler method, although Backward 
Euler method is a stable one for any At. The trapezoidal method 
advantage is that the numerical oscillations appear at the same 
time as the loss in precision, yet the result does not diverge. 
This means that the appearance of the numerical oscillations 
provides a warning that the solution is incorrect, so that the 
calculation can be repeated with a smaller time step. In effect, 
if there are no numerical oscillations, then the precision is 
satisfactory . The Backward Euler does not give such a simple 
monitor of accuracy. 

7.3. Equivalent Model of Nonlinear Storage Elements 


An element whose charge does not change linearly with the voltage, 
i.e. Q CV, but Q = f (V) , is a nonlinear storage element. The 
calculation of the charge increment AQ, which is a consequence of 
the voltage growth Av in time At, is similar to Newton -Raphson's 
iteration procedure for nonlinear elements. But, instead of 
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Newton-Raphson's iteration which uses the tangent, "regula 
falsi" method is used which employs the secant. In Fig. 29 the 
function Q = f (V) is illustrated graphically. 



Fig. 29. The "regula falsi" method application on 
a nonlinear element Q = f( V). 


The nonlinear element, Q = f (V) , is linearized about the point 

V v n 11 


+ ( 3V>n <Cl - V n> - 


(55) 


The charge growth, AQ = Q n+1 ” Q n , is equal to 


Qi°i “ 0 ,. - dn+1 + I n )At/2 * 


‘n+1 ~n 


(56) 


1 ) 


Charge Q 


Eq. 

for Q 


(58) 

( 1 ) 

n+1 


is given accordinq to Eq 
is obtained applying Q° +1 - 

f G tC • 


(55 ) } Q° +1 used in 

f (V° . , ) . The same holds 
n+l 
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Inserting term (55) into (56) gives 


On (V n + l - V * «S+1 + V 4 ^ 2 ' 


J n+1 - 2/4t <3§>n v n+l " [ 2 ' 4t <§>„ V n + 1 J • < 57 > 


dQ, 


Voltage V° +1 and current I° +1 are the values obtained in the 

first iteration. The second step gives the following term for 

o (1) . 

g n+1 * 
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( 1 ) = 
n+1 


Q n + -2±i 

n „o 


Q^j.1 - Q 


v~ - v 

n+1 n 


- (v i + i - v and 


Q n« - Q n * ttj+l + V 4t/2 ' 


which gives 


n+1 Tt 


2 Cl 


- Q 


V u - V 
n+1 n 


» v 1 

n+1 


■ 2 Cl 

Tt 


- Q 


n 


- V 

n+1 n 


V + I_ 
n n 


(58) 


The iterations are continued until the difference, 

AV = - V® +1 , is within the error tolerance, which, as 

earlier mentioned is usually 50 yV. When calculating complex 
networks, the iterations are done together with Newton-Raphson ' s 
iterations for the nonlinear resistance elements. An equivalent 
circuit for the nonlinear storage element Q = f (V) is given in 
Fig. 30 with the conductance 
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Fig, SO. An equivalent circuit of the nonlinear 
storage element s Q = f (V) . 


Q n+1 " Q n 

The secant slope, ■ — , can be calculated with second 

V° - V 
n+1 n 

degree accuracy using the charge derivatives in places v° +1 

and v s 
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n+1 


- Q 
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n+1 


- V 


1 

1 


/dQ» o 


n+1 


/dQx 

WnJ 


n 


( 59 ) 


In the case of semiconductor devices, capacitances which are, by 
the charge change nature in semiconductors, may be defined as 


C 



( 60 ) 


and the capacitance, C, is known in the form C = f (V) . As in the 
previous case, charge Q is obtained by defining the function 
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Q = f (V) using 


Q = | C dV = | f c (V) dV . (61) 

With the calculated charge Q, the method is continued as described 
in terms (57) and (58) . In this case, when the capacitance is 
already given as the function of voltage, it is convenient to use 
term (59) for the secant slope. Then 


Q 


n+1 


- Q 


n 


V' 


n+1 


- V 


7 (c n+l + C n> 


n 


(62) 


Thus the calculation is simplified and there is no need to 
calculate the integral from term (61) . Terms (57) and (58)' 
become: 
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In the case of a nonlinear inductance element, instead of 

<j> = LI we have 4> = f(I) and the procedure is the same as that 

of nonlinear capacitance. According to this 


<J> - f (I) , 


(o) „ 


n+1 


— <f> + 

Y n 


<al>„ « l0 n 


n+1 


- V 
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Expanding A<f> gives 


(o) _ 


♦« = + V_) At/2 . 


n+1 T n n+1 T v n 


By eliminating the A$> from upper terms we obtain; 


<§$>n <Cl ■ I n ) “ <Cl + V n ,tt/2 ' respectively 
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At 
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The second iteration step r using yields: 


^ (1) = 
n+1 
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( 66 ) 


An equivalent circuit, according to term (66) consists of the 
conductance 
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At 


2( <t>n +l ~ * n ) 

X n+1 - X n 


and the current source I = G_v + I . 

L L n n 


7.4 . Transient nonlinear networks analysis procedure 

This paragraph demonstrates the procedure for the numerical 
calculation of a network's transient analysis according to the 
flow diagram of Fig. 31. The network contains linear and nonlinear 
resistive elements, and linear conductances and inductances. 

From Fig. 31 the following is apparent: 

a) The network is defined and every nonlinear resistive 
element is replaced with the linear equivalent circuit 
(conductance G and current source I) . Every capacitance 

C, and inductance L, is replaced with an equivalent circuit 
(conductance G and current source I) . 

b) To start a calculation, the initial capacitors' voltages, 
the initial currents through inductances, and the nonlinear 
elements assumed voltages are given. 

c) The nonlinear elements voltage is calculated using Newton- 
Raphson's iterations since in step b) only the assumed 
voltages were given. The initial capacitors' voltages 

and the initial currents through inductances remain 
unchanged. 

d) Next, the time is advanced by At. This is accomplished 
using numerical integration by calculating the voltage 
increment AV c on capacitors and current increment AI L 
through the inductances. After this is done we turn 
back to step c) and calculate new voltage values on the 
resistance elements. 
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The procedure is continued by repeating steps c) and d) until 
the network's time response is calculated for the entire time 
interval. 




Fig . SI. Numerical calculation procedure of nonlinear 
electrical network transient analysis . 
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